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We study the Schlogl model in a vesicle with semi-permeable membrane. The diffusion constant 
takes a smaller value in the membrane region, which prevents the outflow of self-catalytic product. 
A nonequilibrium state is stably maintained inside of the vesicle. Nutrients are absorbed and waste 
materials are exhausted through the membrane by diffusion. It is interpreted as a model of primitive 
metabolism in a cell. 



I. INTRODUCTION 

Pattern formation in reaction-diffusion systems is actively studied as a theme of nonequilibrium sciences Q . Turing 
type patterns @ , target patterns and spiral patterns Q are found in several reaction diffusion systems. More compli- 
cated patterns such as localized patterns Q and the splitting 0, Q were also found experimentally and theoretically. 
The reaction-diffusion systems have become an important model for various nonlinear processes in biological systems. 
Most reaction-diffusion systems are set up to be as homogeneous as possible to simplify the phenomena and the 
analyses. Recently reaction-diffusion systems in more complicated media such as membranes and microemulsions 
have been studied 0, • 

On the other hand, Oparin considered "corecervate" as an origin of the cell, which played an important role of 
prebiotic chemical evolution Q . He pointed out the importance of chemical reactions confined in the cellular structure. 
Dyson pointed out the importance of the jump process to a highly activated state in a certain bistable system as 
the origin of life |l(J . Recently, Noireaux and Libchaber study a cell-like bioreactor using a phospholipid vesicle [ll| . 
These researchers consider that chemical reactions in a confined cell-like structure are an important step to life. 

A cell is separated from the outside region by a membrane, and complicated chemical reactions occur only inside the 
cell. Some materials and reactants are transported through the membrane. The consumption of nutrient and doing 
away of waste matter through the membrane is a basic process of metabolism, and then highly active nonequilibrium 
state is maintained inside of the cell. Inside and outside of the cell, materials are in the aqueous solution, but the 
membrane is composed of lipid. Therefore, the transport coefficients are different at the membrane. To understand 
the confined or localized nonequilibrium states accompanying the transport of materials through the membrane 
qualitatively, we will study a simple reaction-diffusion system called the Schlogl model in a vesicle with a semi- 
permeable membrane. We do not consider complicated reactions including DNA, RNA and so on, and do not consider 
active transport through the membrane, so it might be too simple as a model of a living cell. That is, our model is 
not a realistic model, but a toy model to consider "metabolism" from a view point of dissipative structure far from 
equilibrium. 

II. LOCALIZED NONEQUILIBRIUM STATE IN THE SCHLOGL MODEL 

The Schlogl model is a set of chemical reactions including three kinds of chemicals U,V and W |l2| represented as 

V + 2U 3U, U ^ W. 

We interpret the three chemicals V,U and W respectively as nutrient material, self-catalytic product and waste 
material. The reaction-diffusion equations for u = [U),v = [V] and w = [W] are expressed as 

— = kivu 2 — k 2 u 3 — k^u + k^w + D u V 2 u, 
ot 

dv o , o 

— = -k x vu 2 + k 2 u A + A,VV 
ot 

— — = k^u — k^w + D w V 2 w. (1) 
ot 

where fci,fc 2 ,/c3 a.ndk 4 are rate constants respectively for the reactions V+2U— >3U, 3U— >V+2U, U— >W and W^ 
U, and D u , D V1 D w are diffusion constants for the materials. In usual reaction-diffusion equations, the diffusion 
constants are assumed to be uniform in space. However, we assume a spherical vesicle surrounded by a membrane as 
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FIG. 1: (a) Localized solution u(r),v(r) and w(r) to eq. (3) for Vq = 1, ro = 2, d = 1 and D^o = 0, D„o = -D«>o = 0.05. (b) 
Critical value Dq c for the localized state as a function of Vq. (c) Extended solution u(r),v(r) and w(r) for Do = 0.05 and 

V = 1.3. 



a model of a primitive cell. The diffusion constants inside of the vesicle are different from the bulk region. We assume 
non-uniform diffusion constants as 

D u (r) = D ul , for r < r , r > r + d, 

= D uQ for r < r < r + 

D v (r) = D vl , for r < r , r > r + d, 

= D v0 for r < r < r + d, 

D w (r) = D w i, for r < r 0l r > r + d, 

= D w0 for r < r < r + d, (2) 

where r is the distance from the center of a spherical cell. The value tq + d denotes the radius of the cell, and d 
is the thickness of the membrane. The diffusion constants are uniform inside and outside of the cell. The diffusion 
constants inside of the cell might be different those outside of the cell, but we assume the same values for the sake 
of simplicity. In all simulations in this paper, we have assumed D„i = D v \ = D w i = 1. The diffusion constants are 
assumed to take smaller values in the membrane region ro < r < ro + d. That is, the transport of materials by the 
diffusion is assumed to be harder in the membrane region. 

If the rate constants k\ ~ fc 4 are assumed to be all 1 for simplicity, and the spherical symmetry is further assumed, 
the model equation is simplified as 

du 2 r, 1 d ( _ o du 

— = vu -u -u + w+ -j— D u r — 
ot r z or \ or 

dv 2 , ..3 , 1 9 fn J* dv 



- -vu ~\~ u + -^t- D v r 
at r z or \ Or 

dw Id/ 2 dw \ 

= u- w + ^Q- r { D ™ r gf)- (3) 

The diffusion constants are given by eq. (2). There is a thermal equilibrium state in eq. (3), because the chemical 
reactions are reversible. It is an important point different from the usual reaction diffusion equation, in which some 
irreversible reactions are assumed in most cases. The equilibrium state is characterized by the condition u = v = w. 
In our simulation, we have imposed fixed boundary conditions of v = vq and w = at a radius r = R = 20 outside of 
the cell, and no-flux boundary conditions of du/dr = at r = R = 20, and du/dr — dv/dr = dw/dr — at r = 0. 
The radius ro and the thickness d are fixed to be ro = 2 and d = 1 for most cases. 

The fixed boundary condition v ~ vq at r = R generates a nonequilibrium state. The completely equilibrium 
state satisfies u — v — w, and it is attained only for the case of i>o = wo = 0. If vq is not zero, the system is 
deviated from the equilibrium state. There is another uniform solution u — w — and v — vq. The eigenvalues of 
this uniform state for the uniform perturbation with wavenumber k = is -2,0 and 0. It implies that the uniform 
solution is marginally stable. Another localized solution can appear. Firstly, we show numerical results in case of 
D u i — D v i — D w i — 1, D u q = and D v q = D w q = Dq. Here, we consider a situation that the self-catalytic product 
U cannot diffuse into the membrane and therefore D u q = 0, assuming that U is a material such as an enzyme of 
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FIG. 2: (a) Stationary solution u(r),v(r) and w(r) to eq. (3) for vo = 1, To = 2, d = 1 and D„o = 0.01, D„o = 0, D^o = 0.01 
and Dui = D v i = D w i = 1. (b) Stationary solution u(r),v(r) and w(r) for D u o = 0, D v o = 0.01, L> m o = and D u i = -Dui = 
D w i — 1. (c) Stationary solution w(r),«(r) and w(r) for D u o = O„o = D w o = and D„i = D v i = D m i = 1. The initial values 
are u(r) = 0.6, v(r) = 0.5, w(r) = 0.4 for r < ro + d and u(r) — v(r) = 0, w(r) = wo = 1 for r > ro + rf. 
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FIG. 3: (a) Concentration difference v(0) — w(0) as a function of the cell size ro. The upper branch denotes the localized 
nonequilibrium state and the lower branch denotes an extended nearly-equilibrium solution, (b) Extended nearly-equilibrium 
solution at D u q = 0, D v q — D w q ~ 0.05, Vo ~ 1 and ro = 2. 



protein and the diffusion hardly occurs through the membrane. Numerical simulation was performed with the Eulcr 
method with mesh sizes At = 0.001 and Ar = 0.1. Figure 1(a) is a stationary state for D u q = 0, D v $ = D w q = 0.05 
and vq — 1. The nutrient V is absorbed through the membrane r < r < r + d, the concentration of the self-catalytic 
product U is highly maintained inside of the cell, and the waste material W diffuses out from the membrane and 
it is absorbed at r = R. The equilibrium condition u = v = w is not satisfied inside of the cell. Nonequilibrium 
localized state is maintained inside of the cell, owing to the inward flow of the nutrient and the outward flow of the 
waste material. This can be interpreted as a simple and primitive model of metabolism. The self-catalytic product 
U is not transported by diffusion through the membrane because of D u q — 0, but it is produced by V and therefore 
u changes continuously in space. As the concentration of the nutrient vq is decreased for Do = 0.05, the localized 
state disappears and it decays to the uniform death state u = w = for vq < 0.71. As Dq is increased for vq = 1, the 
localized state disappears for Dq > 0.121 and leads to the death state u = w = 0. Figure 1(b) displays the critical 
value Do c as a function of vq. The localized state is stable below the critical curve. The localized state cannot be 
maintained, if the diffusion constants are uniform D v q — D v \ and D w q — D w \. For vq < 1.1, Dq c increases and 
the localized state is more stable for larger vq. However, for vo > 1.1, Dq c rapidly decreases. In this region, the 
localized state does not decay to the zero state u — but leads to an extended state for Dq > D$ c . An example of the 
extended state is shown in Fig. 1(c) at vq = 1.3 and Do — 0.05. Inside of the cell r < ro + d, the equilibrium condition 
u = v — w is almost satisfied. We have found only the three types of solutions, i.e., the localized state, the uniform 
death state and the extended state, in our numerical simulation. Our localized nonequilibrium state is maintained 
by the outer boundary conditions. Blanchedeau et al. found spatial bistability of the thermodynamic state and the 
flow state in the chlorine-dioxide-iodide reaction in a thin film of gel fed from one side [l3l |. The phenomenon of the 
spatial bistability is also maintained by the supply of chemicals at the boundary. 
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FIG. 4: (a) Localized solution u(r),v(r) and w(r) to eq. (3) for Vo = 1, T"o = 2, d = 1 and D u o = -Duo = = 0.02. (b) 
Critical value Do c for the localized state as a function of vo- (c) Localized solution of u(r) to eq. (5) at a = 0.979, -Do = 0.02 
and «o = 1, and w = au and v = «o — u — w are also plotted. 



If we consider nonzero diffusivity at least for V and W, the reactants and the waste materials are transported 
through the membrane and the nonequilibrium state is maintained. If D v a = and D u o ^ 0, D w o ^ 0, the nutrient 
is not supplied through the membrane, and only the death state u = w — appears as shown in Fig. 2(a). If 
D u q = D w q = and D v q ^ 0, the nutrient is supplied but the product and the waste material are not sent out from 
the cell. Then, an equilibrium state u = v = w = vq is attained for < r < r$ as shown in Fig. 2(b). If the diffusion 
constants D u0 ,D v0 and D w0 are all 0, the cell is completely isolated from the outside. In that case, the equilibrium 
state u = v = w is attained inside of the cell as shown in Fig. 2(c). In this case, the final concentration of u,v and 
w depends on the initial concentration inside of the cell. In the simulation of Fig. 2(c), the initial concentrations 
were u(r) = 0.6, v(r) = 0.5 and w(r) = 0.4 inside of the cell r < ro + d. As a result of the conservation of the total 
concentration, the final concentration becomes the average value 0.5 of the initial concentrations. Thus, we see that 
the localized state shown in Fig. 1(a) is maintained only in the nonequilibrium condition. 

Even if D v q — D w q = Dq << 1 is not zero, the localized nonequilibrium state such as shown in Fig. 1(a) cannot 
be maintained, if the cell size ro + d is sufficiently large, because the central region of the cell is far away from the 
outside and an nearly equilibrium state is realized. Figure 3(a) displays the difference v(0) — w(0) of the concentration 
as a function of ro for Dq = 0.05 and vq = 1. The parameter ro is gradually increased from 1 to 4 for numerical 
simulations shown by rhombi. Then, ro is gradually decreased from 4 to 1, and the results of numerical simulations 
are denoted by crosses. If ro is larger than 3.2, the localized nonequilibrium state disappears abruptly and it leads to 
an extended state which is similar to the state shown in Fig. 3(a). In the extended state, the equilibrium condition 
is nearly satisfied inside of the cell and therefore v(0) — w(0) takes a very small value. The extended state is stable 
even for ro < 3.2, that is, the localized nonequilibrium state and the extended nearly equilibrium state are bistable 
as shown Fig. 3(a). The extended nearly equilibrium state at vo = 1, ro = 2 and Dq = 0.05 is shown in Fig. 3(b), 
which is bistable with the localized state as shown in Fig. 1(a) at the same parameters. Our numerical result is that 
the localized nonequilibrium state exists only for small cell size ro + d. The result suggests that the cell size cannot 
become too large in the primitive type metabolism, because the material transport occurs only through the simple 
diffusion. More elaborate active transport mechanisms such as a blood vessel network are necessary to maintain a 
large body. 

The localized state appears even if D u q is not zero. We consider hereafter the case of D u = D v = D w for the sake 
of simplicity. We assume D u \ = D v \ = D v \ = 1, and D u0 = D v0 = D w0 is expressed as D - Figure 4(a) displays a 
localized state for vq = 1 and Do = 0.02. Nonequilibrium active state is maintained inside of the cell, owing to the 
inward flow of the nutrient and the outward flow of the waste material. Figure 4(b) displays the critical value Dq c as 
a function of vq. The localized state is stable below the critical curve. The phase diagram is qualitatively the same 
as Fig. 1(b). The stable region of the localized state becomes smaller, because the self-catalytic product U can diffuse 
out. When D u {r) = D v (r) = D w (r) — D(r) is satisfied, the sum s — u + v + w obeys 

^ = iA(V^y (4) 

dt r 2 dr \ dr ) 

The sum s therefore approaches a uniform state s(r) = ^o, where i>o is determined by the boundary condition. We 
have checked that the sum s = u + v + w~vois actually satisfied in the stationary state shown in Fig. 4(a). The 
concentration w is almost equal to u as is seen in Fig. 4(a), but w is slightly smaller than u. If w = au (a ~ 1) is 
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assumed, v — vq — (1 + a)u is obtained. Substitution of these relations into eq. (3) yields 



du L d ( Qua \ 



(5) 



Figure 4(c) displays a localized solution of u to eq. (5) for D Q = 0.02, v Q = 1, w — au and v = 1 — u — au. The 
parameter a — 0.979 is assumed, which is numerically obtained from w(0)/u(0) in Fig. 4(a). The profiles of u,v and 
w in Fig. 4(c) are close to those shown in Fig. 4(a). It implies that the assumption w — au is not exactly satisfied in 
the localized solution of eq. (3), but the approximation is rather good. 

III. LOCALIZED SOLUTION IN THE GINZBURG-LANDAU EQUATION WITH INHOMOGENEOUS 

DIFFUSION CONSTANT 

To understand the mechanism of the stability of the localized state, we study the simpler model (5) further. There 
are three uniform stationary solutions u = and u = uq± = (vq — 1 — a ± \J (vq — 1 — a) 2 — 4(1 — a))/2 to eq. (5), if 
(vq — 1 — a) 2 > 4(1 — a). If u/u + is rewritten as a rescaled variable u, the corresponding reaction-diffusion equation 
becomes the Ginzburg-Landau type equation of the form: 



where a = uo-/uo+ and D(r) = 1 for r > ro + d and r < ro and D{r) = Dq for ro < r < ro + d. There are two stable 
uniform solutions u — and u = 1 for < a < 1 in this model. There is a potential function for this model: 



Equation (6) is a reaction diffusion equation with inhomogeneous diffusion constant. Several authors studied pinning 
phenomena of wave propa gat ion in reaction diffusion systems with some inhomogeneities |l4l 1 1 5| and discrete type 
reaction diffusion systems |l6t Il7|. A localized solution to eq. (6) is related to these pinning phenomena. 

When the parameter a — 0.5, the potential takes the same value for the two uniform solutions u = and u = 1. 
Figure 5(a) is a localized state for a = 0.5, r = 2, d = 1 and D = 0.015. As D is increased, the localized state 
disappears and it decays to the zero state u — 0. We have investigated the critical value of Do above which the 
localized state disappears. Figure 5(b) displays the critical value Dq c as a function of the thickness d for a — 0.5 and 
ro + d/2 = 5. The localized state is most stable at d ~ 1.5, which is the optimum thickness of the membrane. It is 
natural that localized activity in a cell with a thinner membrane is unstable, but it is not so trivial that the critical 
Dq c decreases for large d. For a < 0.5, and the potential effect makes the volume of u = 1 increase, because the 
potential is lower for u = 1. If the radius ro is small, the curvature effect tends to make the inner volume of u = 1 
shrink. The localized solution is maintained by the pinning effect. The stability of the localized state is determined 
by the three effects. Figure 5(c) displays the value of 1 — u at r = ro + d as a function of ro for a = 0.1, d = 1 and 
Dq = 0.05. The localized solution disappears at ro = 3, when the radius ro is gradually increased from 1, and changed 
into the uniform state u = 1. It is because the curvature effect becomes weak for larger ro. The uniform state u = 1 
remains stable even for ro < 3, when ro is gradually decreased from ro = 4. This behavior is somewhat similar to 
Fig. 3(a). Figure 5(d) displays the critical value of Do as a function of ro for a = 0.1 and d = 1. Under the critical 
curve, the localized solution is stable. For ro < 1.9, the localized solution changes into the uniform solution u — 0, 
when Dq is larger than Dq c - For ro > 1.95, the localized solution changes into the uniform solution u = 1, when D 
is larger than Dq c . This is because the curvature effect becomes too large for ro < 1.9 and the inner volume of u = 1 
shrinks for Do > Do c , even though the state of u = 1 is preferable in the potential effect. 

If ro is infinity, the problem is equivalent to a one-dimensional model. Then, we study a simpler one-dimensional 
Ginzburg-Landau equation 



For the one-dimensional model, stationary solutions u(x) can be explicitly solved, because u(x) satisfies the equations: 




(6) 





(7) 



fu 
dx 2 



u(u — a)(u — 1), for x < ro, x > tq + d, 




u{u — a){u — 1), for tq < x < rg + d. 



(8) 
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FIG. 5: (a) Stationary solution u(r) to eq. (6) for a — 0.5, ro = 2, d = 1 and Do = 0.015. (b) Critical value Dq c for the 
localized state as a function of d at a = 0.5 and ro + d/2 = 5. (c) The value 1 — u at the cell surface r = ro + d as a function 
of ro for a = 0.1, d = 1 and Do = 0.05. (d) Critical value Dq c for the localized state as a function of ro at a — 0.1 and d = 1 



The solutions to eq. (8) can be expressed using the elliptic functions. The integral of eq. (8) with respect to x yields 

-y/2{E! + tt 4 /4 - (1 + a)u 3 /3 + au 2 /2}, for x < r , 



da 

dx 
du 

dx 
du 

dx 



= - v /{2/D ){E 2 + u 4 /4 - (1 + a)w 3 /3 + au 2 /2}, for r < x < r + d, 



- y/2{E 3 + u 4 /4 - (1 + a)u 3 /3 + au 2 /2}, for i > r + d, 



(9) 



The constants E\ and £3 of the integral are determined as E\ = —(2a — 1)/12 and E% = 0, because the solution 
approaches u — 1 for a; — > — oo and u = for a; — > oo. The constant E 2 must satisfy 



E 2 = {u 4 /4 - (1 + o)«73 + avr/2}(l/D - 1) + £?i/i?o 



(10) 



at x = r , because of the continuity of D(x)(du(x) / dx) at a; = r . The constant E 2 also needs to satisfy 

E 3 = = {u 4 /4 - (1 + a)" 3 /4 + au 2 /2}(D - 1) + ^2^0, (11) 

at x = r'o + d because of the continuity of D(x)(du(x) / dx) at x = ro + d. We have numerically searched the value of 
u(ro) which satisfies the above conditions (10) and (11). If we can find the value of u(ro) at x = ro, we can construct 
a localized solution by solving eq. (9). Figure 6(a) compares the localized solution (dashed line) obtained by solving 
eq. (9) using it(r ) = 0.610335 with the stationary solution (solid line) obtained by a direct numerical simulation of 
eq. (7) for a = 0.6, d = 1,xq = 9.5 and Do = 0.23. The two curves overlap completely and we cannot distinguish the 
difference. There are discontinuity of the derivative du/ dx at x = ro and ro + d, because of the discontinuity of D(x). 
We have numerically searched the value of u(r ) which satisfies the conditions (10) and (11) for various values of Do- 
Figure 6(b) displays the obtained values of it(ro) as a function of Dq for a — 0.6. There are two solutions, and the 
upper branch represents a stable solution which appears in a direct numerical simulation. The lower one is an unstable 
solution. There is a saddle- node bifurcation at Do = Do c ~ 0.235. There is no localized solution for Do > Do c - The 
critical value Do c decreases with a as shown in Fig. 6(c) for d = 1. The rhombi denote the results by the direct 
numerical simulation and the solid line denotes the results of the numerical analysis for the saddle-node bifurcation. 
The two numerical results coincide well. The critical value Do c for a < 0.5 is equal to Do c for a' = 1 — a > 0.5, owing 
to the symmetry of eq. (7) with respect to the change of variables — 1 — u.. 
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FIG. 6: (a) Stationary solution u(r) for a = 0.6, d = 1 and Do = 0.23 by the direct numerical simulation of eq. (7) (solid line) 
and the analysis by eq. (8) (dashed line) (b) u(ro) as a function of Do for a — 0.6 satisfying eq. (10) and (11). (c) Critical value 
Doc for the localized state as a function of a at d = 1. The dashed line is obtained from eq. (12). 



As a — » 0.5, Dq c approaches 1. In the limit of a = 0.5 and Do — 1, the domain wall is expressed as u = 
[1 — tanh{a(x — xq)}]/2, where a = l/y/8 and xq denotes the position of the domain wall. The potential is evaluated 
in the one-dimension system as 

If a — 0.5 << 1,1 — Do << 1 and xrj is close to ro + d/2, the potential energy E(xo) is approximately given by 

E{xo) ~ E{r + d/2) + ^—^-{x - r - d/2) - ^ ~ [tanh{a(r + d - x )} - tanh{a(r - a? )}]. (12) 

12 o 

The dashed line denotes the critical value of Do, above which a local minimum point disappears in the function 
E(xq). In other words, the domain wall can be trapped and pinned at the local minimum point of E(xq) by the 
inhomogeneity of D(x), if Do is smaller than the critical value. The dashed line is tangential to the solid curve in the 
limit of a = 0.5, but the two critical curves separate away rather rapidly, probably because the approximation of the 
domain wall by u = [1 — tanh{a(a; — xo)}]/2 becomes worse as a is increased. 



IV. SUMMARY AND DISCUSSION 



We have performed numerical simulations of Schlogl model in a vesicle with semi-permeable membrane as a primitive 
model of metabolism. High activity is maintained inside the cell, because of the low diffusivity at the membrane 
region. Absorption of nutrients and release of waste materials occur naturally as a result of the transport through the 
membrane, which is characteristic of the metabolism. The transport of material is essential for the metabolism, and 
it is characteristic of noncquilibrium state. The nonequilibrium state is maintained by the fixed boundary conditions 
v = vo and w — at r = R. Too large cell cannot maintain the nonequilibrium active state, because a nearly thermal 
equilibrium condition is satisfied deep inside the cell. 

A theoretical analysis was performed for the Ginzburg-Landau equation with inhomogeneous diffusion constant 
to understand the mechanism of the stabilization of the localized state better. The Ginzburg-Landau equation was 
derived from the Schlogl model using an assumption w — au, although the parameter a ~ w/u is not obtained 
theoretically. The mechanism of the localized state in the one dimensional model is interpreted as a kind of front 
pinning by the inhomogeneity. It might be an interesting result that there appears an optimal parameter value for 
the localized state with respect to the membrane width d as shown in Fig. 5(b). The front pinning by inhomogeneity 
and discreteness were also studied by several authors, however, our main problem is the localized nonequilibrium 
state in the origin Schlogl model including a thermal equilibrium state. The front pinning is a mechanism for the 
nonequilibrium localized state, but the problem should be further considered from a view point of nonequilibrium 
physics. For example, a localized state is possible even for a cell with large ro in the Ginzburg-Landau equation 
with inhomogeneous diffusion constant, if the parameter a in eq. (6) is assumed to be constant, but a localized 
nonequilibrium state is impossible in the original Schlogl model for such a large cell. 
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Our model might be too simple as a model for a living cell, because we do not consider complicated and realistic 
chemical reactions and various types of active transport of material through the membrane and inside of the cell. 
However, it might be an instructive model to consider metabolism and life conceptually from a view point of a 
dissipative structure far from equilibrium. It might be possible in the future to make an experiment of a primitive 
metabolism in an artificial cell. 
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